The search for subsurface hydrocarbon deposits typically involves a sequence of data acquisition, analysis, and interpretation procedures. The data acquisition phase involves use of an energy source to generate signals that propagate into the earth and reflect from various subsurface geologic structures. The reflected signals are recorded by a multitude of receivers on or near the surface of the earth, or in an overlying body of water. The received signals, which are often referred to as seismic traces, consist of amplitudes of acoustic energy which vary as a function of time, receiver position, and source position and, most importantly, vary as a function of the physical properties of the structures from which the signals reflect. The data analyst uses these traces along with a geophysical model to develop an image of the subsurface geologic structures.
The analysis phase involves procedures that vary depending on the nature of the geological structure being investigated, and on the characteristics of the dataset itself. In general, however, the purpose of a typical seismic data processing effort is to produce an image of the geologic structure from the recorded data. That image is developed using theoretical and empirical models of the manner in which the signals are transmitted into the earth, attenuated by the subsurface strata, and reflected from the geologic structures. The quality of the final product of the data processing sequence is heavily dependent on the accuracy of these analysis procedures.
The final phase is the interpretation of the analytic results. Specifically, the interpreter's task is to assess the extent to which subsurface hydrocarbon deposits are present, thereby aiding such decisions as whether additional exploratory drilling is warranted or what an optimum hydrocarbon recovery scenario may be. In that assessment, the interpretation of the image involves a variety of different efforts. For example, the interpreter often studies the imaged results to obtain an understanding of the regional subsurface geology. This may involve marking main structural features, such as faults, synclines and anticlines. Thereafter, a preliminary contouring of horizons may be performed. A subsequent step of continuously tracking horizons across the various vertical sections, with correlations of the interpreted faults, may also occur. As is clearly understood in the art, the quality and accuracy of the results of the data analysis step of the seismic sequence have a significant impact on the accuracy and usefulness of the results of this interpretation phase.
In principle, the seismic image can be developed using a three-dimensional geophysical model of seismic wave propagation, thereby facilitating accurate depth and azimuthal scaling of all reflections in the data. Accurately specified reflections greatly simplify data interpretation, since the interpretational focus can be on the nature of the geologic structure involved and not on the accuracy of the image. Unfortunately, three dimensional geophysical models frequently require intolerably long computation times, and seismic analysts are forced to simplify the data processing effort as much as possible to reduce the burdens of both analysis time and cost.
In addition to the 3-D computation challenge, the analyst faces a processing volume challenge. For example, a typical data acquisition exercise may involve hundreds to hundreds of thousands of source locations, with each source location having hundreds of receiver locations. Because each source-receiver pair may make a valuable contribution to the desired output image, the data handling load (i.e., the input/output data transfer demand) can be a burden in itself, independent of the computation burden.
FIG. 1 depicts a perspective view of a region 20 of the earth for which a geophysical image is desired. On the surface 18 of the earth are shown a number of shot lines 2 along which the seismic data are acquired. As shown in FIG. 1A, shot lines 2 consist of a sequence of positions at which a seismic source 3 is placed and from which seismic signals 5 are transmitted into the earth. Receivers 4 placed along each line receive the signals from each source position after reflection from various subsurface reflectors 6.
A first method of managing the seismic data burdens discussed above involves careful definition of the region over which the data are acquired. Specifically, use of any available preliminary geologic and geophysical information may facilitate the minimization of the surface area over which seismic data may need to be acquired. Such a minimization will directly reduce the amount of data that is ultimately acquired. Furthermore, similarly careful planning of the spacing between shot lines will optimize the analysis effort by reducing data volume. Finally, optimization of the number of sources and receivers that are used, and of the spacing between adjacent source and receiver positions, will also benefit the data analyst.
None of these efforts can be accomplished without a penalty. For example, relatively wide spacing between shot lines, or between sources and receivers, reduce the resolution of the computed seismic image, thus making interpretation more difficult. In addition, complex geologic features may not be resolvable without relatively close spacing. And finally, certain data acquisition exercises, such as in relatively unexplored areas, do not allow optimization of the surface area over which data is to be acquired. As a result, the data handling burden cannot be entirely eliminated through data acquisition planning.
The recent availability of massively parallel processors offers a significant opportunity to seismic data analysts. Massively parallel processors (MPPs) have multiple central processing units (CPUs) which can perform simultaneous computations. By efficient use of these CPUs, the weeks or months previously required for complex analyses can be reduced to a few days, or perhaps a few hours. However, this significant advantage can only be realized if efficient computational algorithms are encoded in the MPP software. Thus, the opportunity MPPs offer seismic data analysts also creates a challenge for the development of suitable computational algorithms that take advantage of the multiple CPUs.
This challenge can be easily discussed by considering the manner in which computational algorithms have most commonly been written for existing seismic analysis routines. Until recently, computers relied on a mode of operation referred to as sequential computing. Sequential computing involves use of analytic routines that perform only a single procedure, or perhaps focus on a single subset of the data or image, at any given time. This is a direct result of a computer having only one CPU. For that reason, the only optimization procedures that can be employed on single CPU computers are those which increase the efficiency of the processing as to the procedure or subset. Because all calculations must ultimately be performed by that single CPU, however, the options for obtaining high performance are innately limited.
On the other hand, the multiple CPU capability of MPPs offers an obvious simultaneous computation advantage. This advantage is that the total time required to solve a computational problem can be reduced by subdividing the work to be done among the various CPUs, provided that the subdivision allows each CPU to perform useful work while the other CPUs are also performing work. Unfortunately, the disadvantage of multiple CPU hardware is that the sequential processing methods that have long been used in software development must be replaced by more appropriate parallelized computing methods. Simply stated, MPPs require that processing methods be developed which make efficient use of the multiple CPU hardware. Ideally, these methods should organize the distribution of work relatively evenly among the processors, and ensure that all processors are performing necessary computations all of the time, rather than awaiting intermediate results from other processors.
The challenge of defining parallelized processing methods, and of optimizing those parallelized methods once defined, is particularly acute in the seismic data processing arena. Seismic data consists generally of a large number of individual traces, each recorded somewhat independently of the other traces. Logically enough, sequential computing methods that require the analytic focus to be placed on a single calculation at a time adapt well to analysis of these independent traces. This is true even though computational bottlenecks may exist. For example, portions of the analytic sequence may require relatively more computation time than other portions, must be completed before other calculations may proceed, or may rely on similar input data as other traces, for example traveltimes. Since no simultaneous computations occur in sequential processing, none of these bottlenecks lead to a reduction in computational efficiency with a single CPU, except as to the total processing time that is required. Except as to that total time requirement, the existence of such computational bottlenecks does not otherwise pose problems for the analyst. To take full advantage of MPP computing capabilities, however, where the goal is to perform simultaneous processing in all CPUs, methods for optimizing the seismic analysis phase by eliminating such bottlenecks must be developed.
This advantage of an MPP becomes clear by considering the limitation which calculation time places on image region size in single CPU computers. Increasing the size of the image, e.g., by expanding the size of cube 20 in FIG. 1, or increasing the amount of data to be processed, e.g., by adding additional sources 3 and receivers 4 to shotlines 2, increase the total computation. That direct impact on calculation time places a heavy burden on seismic analysts to optimize image size, especially since even small image regions may require weeks of computation time on even the highest speed sequential processing computers. In contrast, efficient processing on MPPS, which may have as many as or more than 256 individual CPUs, should only involve minimally lengthened computation times, since each CPU would assume just a fraction, for example 1/256, of the additional work required by the larger region. This potential for scalability of the image region and the work load required in image generation is a principal benefit of MPPs, a benefit that can only be realized if parallelized seismic processing methods allowing such workload scalability are developed.
Basic considerations for determining efficient parallelized seismic processing methods become evident by reconsidering the above review of the seismic analysis process. As noted, the purpose of seismic analysis is to analyze measured seismic data using geophysical models to develop images of the subsurface. Therefore, each of three principal processing components--data, model, and image--may be considered to be a candidate for distributing computational work among the various processors in an MPP. One option for distributing work among the processors would be to assign different groups of the input seismic trace data to different processors. For example, traces may be grouped by source locations, with different processors being assigned different groups. Similarly, the output image could be subdivided and assigned to different processors. Finally, it may also be possible to subdivide the geophysical model used to generate the output image into groupings that can be assigned to the various processors. (That model is generally considered to be embodied in the arithmetic operations required by the mathematical model that is the subject of the processing effort. For example, in seismic analysis the mathematical model is often based on the wave equation). For example, the data may be transformed into the frequency domain, with individual frequencies assigned to individual processors. It may also be possible to develop combinations of these approaches. For example, groups of processors may be assigned collective responsibility for specific frequencies in the model and all depths in the image, while having individual responsibility for specific horizontal locations in the image. The challenge to the seismic data analyst is to determine methods of subdividing the seismic data, model, and image into components that can be assigned to individual processors in the MPP, thus allowing calculations to be performed in each processor independently of other processors. This subdivision of seismic data analysis into individual components is commonly referred to as seismic decomposition.
In a Single Instruction, Multiple Data stream (SIMD) MPP, the processing elements typically perform the same operation on multiple data streams. An example is the CM2, a product of the Thinking Machines Corporation. These kinds of machines typically lack shared memory, i.e., each processor has its own separate memory unit and the information in the memory cannot be directly accessed by other processors. The individual processors typically have limited computing capability and memory. Because of the large number of processing elements and a lack of shared memory, data transfer between the processing elements is a major bottleneck in efficient utilization of the capability of the machines. Even with sophisticated interconnection techniques, such as in a hypercube arrangement, transfer of data between processors is a major factor in the running time of programs.
Other computers have much more powerful elements in arrays of tens or hundreds. The T3D, a product of Cray Research Corporation, is an example of this kind of machine. Besides having individual processing elements that are much more powerful than those in the CM2, the T3D has fewer of the elements and a physically distributed, logically shared memory. This Multiple Instruction Multiple Data stream (MIMD) machine has different elements performing different operations on different parts of the data at the same time. The reduced number of processing elements means that data does not have to be transferred to as many elements as in a SIMD machine. Because of the increased sophistication and cost of the individual elements and because of their fewer numbers, efficient utilization requires that the load on the processing elements be balanced. An additional factor is that each processing element must accommodate a larger subset of the overall data volume; computations that involve sorting of the data could become more complicated.
As would be familiar to those knowledgeable in the art, prestack migration processes are computationally very demanding. As an alternative to prestack migration, a prestack partial migration could be done followed by a poststack time or depth migration. In many cases, this leads to a significant reduction in computation time without a serious loss in accuracy. One common method of prestack partial migration is the dip moveout (DMO) process. The process of DMO correction is discussed at length in Deregowski, S. M. and Rocca, F., "Geometrical Optics and Wave Theory of Constant Offset Sections in Layered Media", 29 Geoph. Prosp. 374 (1981). The DMO operator can be visualized as the result of two steps. The first step is to find the isochron for a given source and receiver position. The term isochron means the set of points such that the traveltime from the source to each of the points and back to the receiver is equal to the observed travel time. Those knowledgeable in the art would recognize that, for a constant velocity medium, the isochron is an ellipsoid with foci at the shot and the receiver. The normal to any point on this ellipsoid defines both a position on the surface of the earth and a zero-offset travel time corresponding to the ray defined by the normal. The DMO operation consists of moving data from an original observation point to a set of such possible surface positions and determining the corresponding zero-offset travel times. For a constant velocity medium, the surface positions fall along the major axis of the traveltime ellipsoid. The Kirchhoff DMO procedure consists of 1) for each input trace the generation of a set of replacement traces at a number of output surface positions followed by 2) the summation, at each output location, of the replacement traces from each input location, resulting in 3) a partially migrated image.
Deregowski and Rocca showed that a velocity independent DMO operator for a common offset section can be expressed as: ##EQU1## where t.sub.0 is zero-offset time, t.sub.n is the normal moveout (NMO) corrected time, h is half the source-receiver offset and x is the migration distance from the midpoint, i.e., the distance between the midpoint of the source-receiver pair and the position of the equivalent zero-offset replacement trace. Zero-offset time is time measured along the replacement trace. The NMO corrected time is obtained by the usual NMO procedure using the root mean square (rms) velocity.
As is known to those familiar with the art, conventional stacking of common midpoint (CDP) gathers that uses a dip corrected rms velocity is unable to deal with reflectors of different dips. The zero offset section obtained by DMO, on the other hand, correctly represents reflections from different dips. Equation 1 defines a 2-D DMO operator that moves energy into a vertical plane through the source and receiver and produces an image that can be processed by post-stack, zero-offset migration methods at a significant saving in computation time.
When the medium does not have constant velocity, the DMO operator is more complicated. A partial solution for the problem of a medium with constant vertical velocity gradient is given in Artley, C., and Hale, D., "Dip-Moveout Processing for Depth Variable Velocity", 59 Geophysics 610 (April 1994). However, Artley and Hale considered only the inline component of the DMO operator. In fact, for the case of a medium with a vertical velocity gradient, the DMO operator, which is generally referred to as the Migration to Zero Offset (MZO) operator, is no longer limited to a vertical plane passing through the shot and receiver positions. A detailed analysis of the problem of MZO for a constant velocity gradient is given by Dietrich, M., and s Cohen, J. K., "Migration to Zero Offset (DMO) for a Constant Vertical Velocity Gradient: An Analytical Formulation", 41 Geophysical Prospecting 621, (1993). The isochron in this case is a quartic polynomial in x, y and z, and the migration to zero offset (MZO) operator is multivalued. One branch of the operator is a saddle-shaped operator related to smaller dips and approximately corresponds to the conventional DMO operator for the case of no velocity gradient. The second branch of the operator is produced from larger dips or overturned reflectors. Both of these branches of the DMO operator have crossline components. Dietrich and Cohen's work has several limitations. First, their operator is based on ray-theoretical considerations. As a result, where the two branches of the operator come together there is an abrupt termination, with rather strong amplitude, of the primary branch of the operator. In addition, even away from the cusp termination, the amplitudes given by Dietrich and Cohen are inconsistent with the amplitudes given by prior art references relating to the problem of amplitudes for the DMO operator.
U.S. Pat. No. 5,285,422 issued to Gonzalez et. al. teaches a method for performing a 3-D DMO in the presence of a vertical velocity variation in the subsurface. In this method, the DMO operator for constant velocity is applied, and followed by the application of a residual operator that has the effect of squeezing the ellipsoid from the first step. However, the method only addresses the primary branch of the multivalued operator. In U.S. Pat. No. 4,974,212, Scheiman deals with the problem of smearing caused by the operation of forcing the output location into a discrete set of bins on an areal grid. However, that method does not address the problem of vertical inhomogeneities. In addition to these limitations both Gonzalez et. al. and Sheiman both are single CPU-oriented methods and cannot be efficiently implemented on an MPP.
The implementation of a 3-D DMO process on a shared memory parallel computer is discussed by Lu, L., Bell, L., and Lim, K., "Parallel Implementation of 3-D DMO on Shared Memory Systems", 64 SEG Extended Abstracts 214, (1994). Lu et. Al. implement a parallelization of tasks over input shot gathers only, but not in the output phase of the DMO calculations. However, the disclosure is limited to DMO of single offset data, and in addition has several limitations that reduce the efficiency of the overall process. Specifically, because multiple processors may be attempting to write (or add) their output to the same memory location of the MPP, a "memory locking" procedure is incorporated in the process. This memory locking procedure results in a slowdown of the process, since processors will often be delayed in accessing the desired memory locations. In addition, after processing of each group of shots is completed, the disk file is updated using an analogous "file locking". This too is a bottleneck in the efficiency of the method.
An MPP implementation for a 3-D DMO process on a CM2 (or on the Thinking Machine CM5 computer) is discussed by Biondi, B. and Moorhead, W. D., "Data Parallel Algorithm for Kirchhoff 3-D DMO", 62 SEG Extended Abstracts 330, (1992). The implementation has parallelism over input shot gathers and also over output swaths. Successive shot gathers are loaded onto consecutive processors until all the processors are full, each processor having one or more gathers as part of its input. Each processor is also responsible for a number of output locations. During the computation phase, each processor first computes the DMO replacement traces for its input traces and only for the locations for which it also has the output task. Once all the processors have completed calculations for the initial input traces, each processor sends its input traces to the next processor and the DMO operations are repeated. This "pipelining" is repeated, as many times as necessary, until output traces have been computed for all output locations. The pipelining is necessary because of the limited memory available on the individual processors of a CM-2 or CM-5. The processors also drop input traces incrementally, beginning with the nearest offset, to reduce the data transfer. This method requires that the shot gathers be contiguous and in physical order, or the load balancing between the processors will be poor. In addition, since each processor is responsible for contiguous portions of the output, bottlenecks may occur in certain portions of the model where, because of complicated velocity models, more computations are needed. Finally, the load balancing in such a scheme will deteriorate when, as is commonly the case, certain traces are missing from the input data set.
Gonzalez et. Al., Lu et. al. and Biondi and Moorhead deal only with the kinematics of the DMO operation, i.e., the DMO calculation involves a straightforward summation of the generated replacement traces. They do not properly handle the amplitude and phase information in the seismic reflections. The problem of preserving amplitudes and phase information as part of the DMO process is discussed by Black, J. L., Schleicher, K. L. and Zhang, L., "True-Amplitude Imaging and Dip Moveout", 58 Geophysics 47 (January 1993). Black et. al. point out that merely summing the contributions, at the appropriate times, from various input traces can give a misleading image of the subsurface. In particular, the amplitudes of steeply dipping reflectors would be attenuated with respect to horizontal reflectors. Those knowledgeable in the art would recognize that amplitude and phase information play an important role in interpretation of subsurface lithology and fluid content and would hence recognize the desirability of preserving amplitude and phase information. Black et. al. also discuss the problem of aliasing wherein, due to too large a distance between input trace positions, steep dips are improperly imaged. Black et. al. do not, however, suggest a solution to that problem.
Thus, there is a need for an invention that implements an MZO process on an MPP which efficiently load balances the DMO calculations, parallelizes on input and output, is able to handle multiple offsets, is able to handle vertical velocity variations, produces accurate images using both branches of the MZO operator, and does not require the data to be contiguous. In addition, the invention should more accurately preserve amplitude and phase information in the seismic data than does the prior art, and should also avoid aliasing errors. The present invention satisfies this need.